rm(list = ls())
source("code/binary-fit.R")
library("EBMAforecast")
library("xtable")
library("magrittr")
library(tidyr)
####################################################################################
## Table 3
####################################################################################
model_names <- c("空间扩散模型", "时间依赖模型", "地理环境模型",
                 "时空基础模型", "时空全模型")

load("results/ensemble.rda")
load("results/pr_calib.rda")
load("results/pr_test.rda")
load("results/pr_test6.rda")
# take out months with window edge in future
pr_test6 <- pr_test6[!is.na(pr_test6$y), ]
# Combine fit statistics for forecast
tab_ebma <- data.frame(
              Model = c(model_names, "Ensemble"),
              W = c(round(ensemble@modelWeights, 2), ""),
              a0 = c(ensemble@modelParams[1, , 1], NA),
              a1 = c(ensemble@modelParams[2, , 1], NA),
              Brier    = apply(pr_calib[, grep("i[0-9]|ebma", colnames(pr_calib))], 2, brier,   obs=pr_calib$y),
              AUC_ROC  = apply(pr_calib[, grep("i[0-9]|ebma", colnames(pr_calib))], 2, auc_roc, obs=pr_calib$y),
              AUC_PR   = apply(pr_calib[, grep("i[0-9]|ebma", colnames(pr_calib))], 2, auc_pr,  obs=pr_calib$y),
              Brier2   = apply(pr_test6[, grep("i[0-9]|ebma", colnames(pr_test6))], 2, brier,   obs=pr_test6$y),
              AUC_ROC2 = apply(pr_test6[, grep("i[0-9]|ebma", colnames(pr_test6))], 2, auc_roc, obs=pr_test6$y),
              AUC_PR2  = apply(pr_test6[, grep("i[0-9]|ebma", colnames(pr_test6))], 2, auc_pr,  obs=pr_test6$y),
              row.names = NULL
  )
ens_row <- match("Ensemble", tab_ebma$Model)
tab_ebma <- rbind(tab_ebma[ens_row, ], tab_ebma[-ens_row, ])
tab_ebma %>% 
      xtable(digits=c(0, 0, 2, 2, 2, 5, 3, 3, 4, 3, 3)) %>% 
      print(include.rownames=FALSE, file = "figs/table3.html",type = "html")
####################################################################################
  
  
  
####################################################################################
### Figure 5
####################################################################################
  
tab_ebma_long <- tab_ebma %>%
             mutate(W = c(NA, ensemble@modelWeights))
tab_ebma_long %<>% gather(fit_stat, value, Brier:AUC_PR2)
tab_ebma_long %<>% filter(Model != "Ensemble")
tab_ebma_long %<>% filter(fit_stat %in% c("Brier", "AUC_ROC", "AUC_PR"))
tab_ebma_long$fit_stat <- gsub("\\_", "-", tab_ebma_long$fit_stat)
tab_ebma_long$fit_stat <- factor(tab_ebma_long$fit_stat, levels = unique(tab_ebma_long$fit_stat))

p <- ggplot(tab_ebma_long, aes(x = W, y = value)) +
        geom_smooth(method="lm", se=FALSE, color="gray80", alpha=0.3) +
        geom_point(colour="gray80", alpha=0.3) +
        geom_text(aes(label=Model), size=3, alpha = 0.9, position="dodge") +
        facet_wrap( ~ fit_stat, scales="free") +
        theme_bw() + 
        scale_x_continuous(limits = c(0, 0.5),expand = c(.1,0.1)) + 
        labs(y = "")+theme(text = element_text(size=10, family = 'KaiTi'))
ggsave(plot=p, file="figs/Figure_5.jpeg", width=6, height=2, scale=1.5, dpi = 400)
####################################################################################

  
  